#This program solves the Navier Stokes equation. It uses the method of
#constructed solutions to test for convergence, which should be 2nd order 

#from __future__ import division
import numpy as np
import numpy.fft as fft
from numpy import pi 
from pyfft.cl import Plan
import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.elementwise import ElementwiseKernel
import time
import sys
from pyvisfile.vtk import write_structured_grid
from pyvisfile.vtk import (
           UnstructuredGrid, DataArray,
           AppendedDataXMLGenerator,
           VTK_VERTEX, VF_LIST_OF_VECTORS, VF_LIST_OF_COMPONENTS)
#np.set_printoptions(threshold=sys.maxint)
from warnings import filterwarnings
filterwarnings("ignore", module="pyopencl.cache", lineno=336)
filterwarnings("ignore", module="pyopencl", lineno=163)
TIMING = 0 
WGPSIZE_1 = 16 
WGPSIZE_2 = 16 
b_WGPSIZE = 4   
PR_INT = 100
SAVE = 1
start_size = 6 
end_size = 6 

def spread_routine(q, firstn1, nextn1, firstn2, nextn2, F_1, F_2, Mu, Rho, X_1, X_2, N, dx, dtheta, f_1, f_2, mu_1, mu_2, rho_1, rho_2):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            spread_kernel(q, (N/4, N/4), (WGPSIZE_1, WGPSIZE_2), firstn1.data, nextn1.data, firstn2.data, nextn2.data, F_1.data, F_2.data, Mu.data, Rho.data, X_1.data, X_2.data, N, dx, dtheta, shift0, shift1, f_1.data, f_2.data, mu_1.data, mu_2.data, rho_1.data, rho_2.data)

def update_routine(q,  firstn1, nextn1, firstn2, nextn2, X_1, X_2, N, dx):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            updaten_kernel(q, (N/4, N/4), (WGPSIZE_1, WGPSIZE_2), firstn1.data, nextn1.data, firstn2.data, nextn2.data, X_1.data, X_2.data, N, dx, shift0, shift1)

def calc_transf(dx):
    N = np.int(np.round(1/dx,0))
    k = np.arange(0, N/2+1)
    k = np.hstack([k, np.arange(-N/2+1, 0)])
    k1 = np.tile(k, (N, 1))
    k2 = k1.T.copy()
    L_tf = -4./dx**2 * (np.sin(pi*k1/N)**2+np.sin(pi*k2/N)**2) 
    L_tf[0,0] = 1. #avoid dividing by zero
    return L_tf

def calc_inds(N): #calculate shifted indices
    a1 = range(0,N)
    i_1 = np.tile(a1, (N, 1))
    i_2 = i_1.T.copy()
    a2 = np.hstack([range(1, N), 0])
    ip_1 = np.tile(a2, (N,1))
    ip_2 = ip_1.T.copy()
    a3 = np.hstack([N-1, range(0, N-1)])
    im_1 = np.tile(a3, (N,1))
    im_2 = im_1.T.copy()
    return i_1, i_2, ip_1, ip_2, im_1, im_2

def solve(t_step):
    global dx, dt, mu_0, rho_0, b_q, dtheta, mu_bar, rho_bar, dr, N_r, body_f, k_teth, ks_teth, kb_teth
    N = np.int(np.round(1/dx,0))
    copy_gpu(p_gpu, p_last_gpu)

    #print 'u_1 is', u_1_gpu.get().max()
    #print 'u_2 is', u_2_gpu.get().max()
    imag_gpu(u_1_gpu, u_1_tf_gpu)
    imag_gpu(u_2_gpu, u_2_tf_gpu)

    plan.execute(u_1_tf_gpu.data)
    plan.execute(u_2_tf_gpu.data)
    
    interp_kernel(queue, (b_q, 1), (b_WGPSIZE, 1), u_1_gpu.data, u_2_gpu.data, X_1_gpu.data, X_2_gpu.data, N, dx, U_1_gpu.data, U_2_gpu.data)
    add_gpu(X_1_gpu, U_1_gpu, dt, X_1_tilde_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt, X_2_tilde_gpu)

    #calc_lagF_gpu(X_1_gpu, X_2_gpu, rr_lag_gpu, ip_s_lag_gpu, im_s_lag_gpu, ds, F_1_gpu, F_2_gpu)
    zero_gpu(F_1_gpu)
    zero_gpu(F_2_gpu)
    calc_lagF_kernel(queue, (N_r, 1), (b_WGPSIZE, 1), X_1_gpu.data, X_2_gpu.data, X_1_teth_gpu.data, X_2_teth_gpu.data, i_teth_gpu.data, i_t_gpu.data, ip_t_gpu.data, im_t_gpu.data, k_teth, ks_teth, kb_teth, dr, F_1_gpu.data, F_2_gpu.data)
    #print F_1_gpu.get()
    #print F_2_gpu.get()
    scale_gpu(F_1_gpu, dr/dtheta, F_1_gpu)
    scale_gpu(F_2_gpu, dr/dtheta, F_2_gpu)

    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)
    zero_gpu(u_1_tilde_gpu)
    zero_gpu(u_2_tilde_gpu)
    zero_gpu(mu_1_gpu)
    zero_gpu(mu_2_gpu)
    zero_gpu(rho_1_gpu)
    zero_gpu(rho_2_gpu)

    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, Mu_gpu, Rho_gpu, X_1_gpu, X_2_gpu, N, dx, dtheta, u_1_tilde_gpu, u_2_tilde_gpu, mu_1_gpu, mu_2_gpu, rho_1_gpu, rho_2_gpu)
    add_gpu(u_1_tilde_gpu, body_gpu, body_f, u_1_tilde_gpu)
    #print mu_1_gpu.get().max()
    #print rho_1_gpu.get().max()
    add_scalar_gpu(mu_1_gpu, -mu_bar, mu_1_gpu)
    add_scalar_gpu(mu_2_gpu, -mu_bar, mu_2_gpu)
    add_scalar_gpu(rho_1_gpu, -rho_bar, rho_1_gpu)
    add_scalar_gpu(rho_2_gpu, -rho_bar, rho_2_gpu)
    #check power identity
    #print 'power1 is', (U_1_gpu.get()*Rho_gpu.get()+U_2_gpu.get()*Rho_gpu.get()).sum()*dtheta-(u_1_gpu.get()*rho_1_gpu.get()+u_2_gpu.get()*rho_2_gpu.get()).sum()*dx*dx
    #print 'power2 is', (U_1_gpu.get()*F_1_gpu.get()+U_2_gpu.get()*F_2_gpu.get()).sum()*dtheta-(u_1_gpu.get()*u_1_tilde_gpu.get()+u_2_gpu.get()*u_2_tilde_gpu.get()).sum()*dx*dx
    #check momentum conservation
    #print 'mom1 is', (u_1_gpu.get()*(rho_1_gpu.get()+rho_0)).sum()*dx*dx
    #print 'mom2 is', (u_2_gpu.get()*(rho_2_gpu.get()+rho_0)).sum()*dx*dx
    #print 'mom1 is', u_1_gpu.get().sum()*dx*dx
    #print 'mom2 is', u_2_gpu.get().sum()*dx*dx
    if PR_INT:
        vort_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, dx, N, vort_gpu.data)

    add_gpu(f_1_gpu, u_1_tilde_gpu, 1, f_1_gpu) 
    add_gpu(f_2_gpu, u_2_tilde_gpu, 1, f_2_gpu) 

    vc_laplace_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, mu_1_gpu.data, mu_2_gpu.data, dx, N, u_1_tilde_gpu.data, u_2_tilde_gpu.data) 

    add_gpu(f_1_gpu, u_1_tilde_gpu, 1.0, f_1_gpu) 
    add_gpu(f_2_gpu, u_2_tilde_gpu, 1.0, f_2_gpu) 

    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, dx, N, u_1_tilde_gpu.data, u_2_tilde_gpu.data)
    #zero_gpu(u_1_tilde_gpu) #turn off convection
    #zero_gpu(u_2_tilde_gpu)

    add_gpu(f_1_gpu, u_1_tilde_gpu, -rho_0, f_1_gpu)
    add_gpu(f_2_gpu, u_2_tilde_gpu, -rho_0, f_2_gpu)

    dal_term_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, u_1_o_gpu.data, u_2_o_gpu.data, u_1_tilde_gpu.data, u_2_tilde_gpu.data, rho_1_gpu.data, rho_2_gpu.data, dt, N, mu_1_gpu.data, mu_2_gpu.data)

    copy_gpu(u_1_gpu, u_1_o_gpu)
    copy_gpu(u_2_gpu, u_2_o_gpu)

    copy_gpu(f_1_gpu, u_1_tilde_gpu)
    copy_gpu(f_2_gpu, u_2_tilde_gpu)

    add_gpu(f_1_gpu, mu_1_gpu, -1.0, f_1_gpu)
    add_gpu(f_2_gpu, mu_2_gpu, -1.0, f_2_gpu)
    ##########################################################solve for u_tilde
    div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), f_1_gpu.data, f_2_gpu.data, dx, N, p_gpu.data)
    imag_gpu(p_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(p_tf_gpu, L_tf_gpu, p_tf_gpu)
    p_filt_kernel(queue, (1,1), (1,1), p_tf_gpu.data)

    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)
    #print 'p is', p_gpu.get().max()

    grad_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), p_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1, f_2_gpu)

    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tf_gpu.data, u_2_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, dt, rho_0, mu_0, N, f_1_tf_gpu.data, f_2_tf_gpu.data)

    nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_1_tf_gpu.data, N)
    nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_2_tf_gpu.data, N)

    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)

    add_gpu(f_1_gpu, u_1_tilde_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_tilde_gpu, 0.5, f_2_gpu)

    plan.execute(f_1_tf_gpu.data, inverse=True)
    plan.execute(f_2_tf_gpu.data, inverse=True)
    
    real_gpu(f_1_tf_gpu, u_1_tilde_gpu)
    real_gpu(f_2_tf_gpu, u_2_tilde_gpu)
    #print 'u_1 is', u_1_tilde_gpu.get().max()
    #print 'u_2 is', u_2_tilde_gpu.get().max()
    ############################################################

    ############################################################solve for u_soln

    #calc_lagF_gpu(X_1_tilde_gpu, X_2_tilde_gpu, rr_lag_gpu, ip_s_lag_gpu, im_s_lag_gpu, ds, F_1_gpu, F_2_gpu)
    zero_gpu(F_1_gpu)
    zero_gpu(F_2_gpu)
    calc_lagF_kernel(queue, (N_r, 1), (b_WGPSIZE, 1), X_1_tilde_gpu.data, X_2_tilde_gpu.data, X_1_teth_gpu.data, X_2_teth_gpu.data, i_teth_gpu.data, i_t_gpu.data, ip_t_gpu.data, im_t_gpu.data, k_teth, ks_teth, kb_teth, dr, F_1_gpu.data, F_2_gpu.data)
    scale_gpu(F_1_gpu, dr/dtheta, F_1_gpu)
    scale_gpu(F_2_gpu, dr/dtheta, F_2_gpu)

    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx)

    zero_gpu(u_1_gpu)
    zero_gpu(u_2_gpu)
    zero_gpu(mu_1_gpu)
    zero_gpu(mu_2_gpu)

    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, Mu_gpu, Rho_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx, dtheta, u_1_gpu, u_2_gpu, mu_1_gpu, mu_2_gpu, rho_1_gpu, rho_2_gpu)
    add_gpu(u_1_gpu, body_gpu, body_f, u_1_gpu)
    add_scalar_gpu(mu_1_gpu, -mu_bar, mu_1_gpu)
    add_scalar_gpu(mu_2_gpu, -mu_bar, mu_2_gpu)
    add_scalar_gpu(rho_1_gpu, -rho_bar, rho_1_gpu)
    add_scalar_gpu(rho_2_gpu, -rho_bar, rho_2_gpu)
    if SAVE:
        copy_gpu(mu_1_gpu+mu_0, mu_1_save_gpu)
        copy_gpu(mu_2_gpu+mu_0, mu_2_save_gpu)

    scale_gpu(rho_1_gpu, 0.5, rho_1_gpu)
    scale_gpu(rho_2_gpu, 0.5, rho_2_gpu)
    if SAVE:
        copy_gpu(rho_1_gpu, rho_1_save_gpu)
        copy_gpu(rho_2_gpu, rho_2_save_gpu)
        add_scalar_gpu(rho_1_save_gpu, rho_0, rho_1_save_gpu)
        add_scalar_gpu(rho_2_save_gpu, rho_0, rho_2_save_gpu)

    add_gpu(f_1_gpu, u_1_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, 0.5, f_2_gpu)

    vc_laplace_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, mu_1_gpu.data, mu_2_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data) 

    add_gpu(f_1_gpu, u_1_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, 0.5, f_2_gpu)

    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data)
    #zero_gpu(u_1_gpu) #turn off convection
    #zero_gpu(u_2_gpu)

    add_gpu(f_1_gpu, u_1_gpu, -rho_0*0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -rho_0*0.5, f_2_gpu)

    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_o_gpu.data, u_2_o_gpu.data, dx, N, mu_1_gpu.data, mu_2_gpu.data)
    #zero_gpu(mu_1_gpu) #turn off convection
    #zero_gpu(mu_2_gpu)

    add_gpu(u_1_gpu, mu_1_gpu, 1.0, u_1_gpu)
    add_gpu(u_2_gpu, mu_2_gpu, 1.0, u_2_gpu)

    scale_gpu(u_1_gpu, 0.5, u_1_gpu)
    scale_gpu(u_2_gpu, 0.5, u_2_gpu)

    dal_term_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, u_1_o_gpu.data, u_2_o_gpu.data, u_1_gpu.data, u_2_gpu.data, rho_1_gpu.data, rho_2_gpu.data, dt, N, u_1_gpu.data, u_2_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1.0, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1.0, f_2_gpu)

    div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), f_1_gpu.data, f_2_gpu.data, dx, N, p_gpu.data)
    imag_gpu(p_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(p_tf_gpu, L_tf_gpu, p_tf_gpu)
    p_filt_kernel(queue, (1,1), (1,1), p_tf_gpu.data)
    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)
    #print 'p is', p_gpu.get().max()

    grad_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), p_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1, f_2_gpu)

    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tf_gpu.data, u_2_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, dt, rho_0, mu_0, N, f_1_tf_gpu.data, f_2_tf_gpu.data)

    nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_1_tf_gpu.data, N)
    nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_2_tf_gpu.data, N)

    plan.execute(f_1_tf_gpu.data, inverse=True)
    plan.execute(f_2_tf_gpu.data, inverse=True)
    
    real_gpu(f_1_tf_gpu, u_1_gpu)
    real_gpu(f_2_tf_gpu, u_2_gpu)

    interp_kernel(queue, (b_q, 1), (b_WGPSIZE, 1), u_1_gpu.data, u_2_gpu.data, X_1_tilde_gpu.data, X_2_tilde_gpu.data, N, dx, F_1_gpu.data, F_2_gpu.data)
    add_gpu(U_1_gpu, F_1_gpu, 1, U_1_gpu)
    add_gpu(U_2_gpu, F_2_gpu, 1, U_2_gpu)
    add_gpu(X_1_gpu, U_1_gpu, dt/2, X_1_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt/2, X_2_gpu)
    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_gpu, X_2_gpu, N, dx)
    ##############################################################
    return

def err(u_1, u_2, u_1_o, u_2_o, X_1, X_2, X_1_o, X_2_o, p, p_o):
    global dx, dtheta
    NN_o = u_1_o.shape[0]
    errinf1 = 0.
    errinf2 = 0.
    err11 = 0.
    err12 = 0.
    err21 = 0.
    err22 = 0.

    NN_r_o = 3*NN_o/8
    NN_s_o = 75*NN_o/16
    NN_r = 3*NN_o/4
    NN_s = 75*NN_o/8
    X_1_o = X_1_o.reshape(NN_s_o,NN_r_o)
    X_2_o = X_2_o.reshape(NN_s_o,NN_r_o)
    X_1 = X_1.reshape(NN_s,NN_r)
    X_2 = X_2.reshape(NN_s,NN_r)
    errinfX = 0.
    err1X = 0.
    err2X = 0.

    errinfp = 0
    err1p = 0
    err2p = 0

    for i in range(0,NN_o):
        for j in range(0,NN_o):
            if (np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))>errinf1):
                errinf1 = np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))
            if (np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))>errinf2):
                errinf2 = np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))
            err21 += np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))**2
            err22 += np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))**2
            err11 += np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))
            err12 += np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))
            if (np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))>errinfp):
                errinfp = np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))
            err2p += np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))**2
            err1p += np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1])) 

    err21 = (dx/2)*np.sqrt(err21)
    err22 = (dx/2)*np.sqrt(err22)
    err11 = (dx**2/4)*err11
    err12 = (dx**2/4)*err12
    err2p = (dx/2)*np.sqrt(err2p)
    err1p = (dx**2/4)*err1p

    for i in range(0,NN_r_o):
        for j in range(0,NN_s_o):
            if (np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)>errinfX):
                errinfX = np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)
            err2X += np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)**2
            err1X += np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)

    err2X = np.sqrt(dtheta/4)*np.sqrt(err2X)
    err1X = (dtheta/4)*err1X

    return errinf1, errinf2, err11, err12, err21, err22, errinfX, err1X, err2X, errinfp, err1p, err2p 

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
Nyq_gpu = ElementwiseKernel(ctx,
	"double2 *u, double *filt, double2 *f",
	"f[i] = filt[i]*u[i]",
	"Nyq_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

calc_lagF_gpu = ElementwiseKernel(ctx,
	"double *X_1, double *X_2, "
	"double *rr_lag, unsigned *ip_s_lag, unsigned *im_s_lag, " 
	"double ds, "
	"double *F_1, double *F_2",
	"F_1[i] = 1./(ds*ds)*(1.+sin(2*pi*rr_lag[i]-pi/2))"
	   "*(X_1[ip_s_lag[i]]-2.*X_1[i]+X_1[im_s_lag[i]]); "
	"F_2[i] = 1./(ds*ds)*(1.+sin(2*pi*rr_lag[i]-pi/2))"
	   "*(X_2[ip_s_lag[i]]-2.*X_2[i]+X_2[im_s_lag[i]])",
	#"F_1[i] = 1./(ds*ds)"
	#   "*(X_1[ip_s_lag[i]]-2.*X_1[i]+X_1[im_s_lag[i]]); "
	#"F_2[i] = 1./(ds*ds)"
	#   "*(X_2[ip_s_lag[i]]-2.*X_2[i]+X_2[im_s_lag[i]])",
        #"F_1[i] = 0.;"
        #"F_2[i] = 0.",
	"calc_lagF_gpu",
	preamble="""
	#define pi (double)3.141592653589793
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

complex_mult_gpu = ElementwiseKernel(ctx,
	"double2 *u_1, double2 *u_2, double2 *f",
	"f[i] = complex_mul(u_1[i], u_2[i])",
	"complex_mult_gpu",
	preamble="""
	#define complex_ctr(x, y) (double2)(x, y)
	#define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

mult_gpu = ElementwiseKernel(ctx,
	"double *u_1, double *u_2, double *f",
	"f[i] = u_1[i]*u_2[i]",
	"mult_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

add_gpu = ElementwiseKernel(ctx,
	"double *u_1, double *u_2, double a, double *f",
	"f[i] = u_1[i]+a*u_2[i]",
	"add_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")
    
add_scalar_gpu = ElementwiseKernel(ctx,
	"double *u, double a, double *f",
	"f[i] = u[i]+a",
	"add_scalar_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")
    
scale_gpu = ElementwiseKernel(ctx,
	"double *u, double a, double *f",
	"f[i] = a*u[i]",
	"scale_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

copy_gpu = ElementwiseKernel(ctx,
	"double *u, double *f",
	"f[i] = u[i]",
	"copy_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

real_gpu = ElementwiseKernel(ctx,
	"double2 *x, double *z",
	"z[i] = x[i].x",
	"real_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

imag_gpu = ElementwiseKernel(ctx,
	"double *x, double2 *z",
	"z[i].x = x[i];"
	"z[i].y = 0",
	"imag_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

zero_gpu = ElementwiseKernel(ctx,
	"double *x",
	"x[i] = 0",
	"zero_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")
 
null_imag_gpu = ElementwiseKernel(ctx,
	"double2 *x, double2 *z",
	"z[i].x = x[i].x;"
	"z[i].y = 0",
	"null_imag_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

calc_p_gpu = ElementwiseKernel(ctx,
	"double2 *w_tf, double2 *L_tf, double2 *p_tf",
	"p_tf[i] = complex_div(w_tf[i], L_tf[i])",
	"calc_p_gpu",
	preamble="""
	#define complex_ctr(x, y) (double2)(x, y)
	#define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

if False:
    solve_gpu = ElementwiseKernel(ctx,
    "double2 *u_1_tf, double2 *u_2_tf, double2 *L_tf,"
    "double2 *f_1_tf, double2 *f_2_tf, double2 *dp_1_tf,"
    "double2 *dp_2_tf, double dt, double rho, double mu,"
    "double2 *u_1_tf_tilde, double2 *u_2_tf_tilde",
    "double2 rhs_1 = complex_mul(1+mu*dt/(2*rho)*L_tf[i], u_1_tf[i])+dt/rho*(f_1_tf[i]-dp_1_tf[i]);"
    "double2 rhs_2 = complex_mul(1+mu*dt/(2*rho)*L_tf[i], u_2_tf[i])+dt/rho*(f_2_tf[i]-dp_2_tf[i]);"
    "double2 lhs = 1-mu*dt/(2*rho)*L_tf[i];"
    "u_1_tf_tilde[i] = complex_div(rhs_1, lhs);"
    "u_2_tf_tilde[i] = complex_div(rhs_2, lhs);",
    "solve_gpu",
    preamble="""
    #define complex_ctr(x, y) (double2)(x, y)
    #define complex_mul(a, b) complex_ctr(mad(-(a).y, (b).y, (a).x * (b).x), mad((a).y, (b).x, (a).x * (b).y))
    #define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
    """)

ib_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    #pragma OPENCL EXTENSION cl_khr_fp64 : enable
    //#pragma OPENCL EXTENSION cl_amd_fp64 : enable
    __kernel void interp(__global const double *u_1, __global const double *u_2, __global const double *X_1, __global const double *X_2, const unsigned N, const double dx, __global double *U_1, __global double *U_2)
    {
      unsigned k = get_global_id(0);
      double s_1_stag = (X_1[k]+0.5)/dx;
      double s_2_stag = (X_2[k]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      double r_1_stag = s_1_stag-i_1_stag;
      double r_2_stag = s_2_stag-i_2_stag;
      double s_1_c = (X_1[k]+0.5-dx/2)/dx;
      double s_2_c = (X_2[k]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      double r_1_c = s_1_c-i_1_c;
      double r_2_c = s_2_c-i_2_c;
      
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;
      ind_1x = (i_1_stag-1 & N-1);
      ind_1y = (i_1_stag & N-1);
      ind_1z = (i_1_stag+1 & N-1);
      ind_1w = (i_1_stag+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      double sb_tests0, sb_tests1, sb_tests2, sb_tests3;
      double sb_tests4, sb_tests5, sb_tests6, sb_tests7;
      double sb_tests8, sb_tests9, sb_testsa, sb_testsb;
      double sb_testsc, sb_testsd, sb_testse, sb_testsf;

      sb_tests0 = u_1[ind_1x+N*ind_2x];
      sb_tests1 = u_1[ind_1y+N*ind_2x];
      sb_tests2 = u_1[ind_1z+N*ind_2x];
      sb_tests3 = u_1[ind_1w+N*ind_2x];
      sb_tests4 = u_1[ind_1x+N*ind_2y];
      sb_tests5 = u_1[ind_1y+N*ind_2y];
      sb_tests6 = u_1[ind_1z+N*ind_2y];
      sb_tests7 = u_1[ind_1w+N*ind_2y];
      sb_tests8 = u_1[ind_1x+N*ind_2z];
      sb_tests9 = u_1[ind_1y+N*ind_2z];
      sb_testsa = u_1[ind_1z+N*ind_2z];
      sb_testsb = u_1[ind_1w+N*ind_2z];
      sb_testsc = u_1[ind_1x+N*ind_2w];
      sb_testsd = u_1[ind_1y+N*ind_2w];
      sb_testse = u_1[ind_1z+N*ind_2w];
      sb_testsf = u_1[ind_1w+N*ind_2w];

      double q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      double q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      double v1x, v1y, v1z, v1w;
      double v2x, v2y, v2z, v2w;

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      double ws0, ws1, ws2, ws3;
      double ws4, ws5, ws6, ws7;
      double ws8, ws9, wsa, wsb;
      double wsc, wsd, wse, wsf;
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      ws0 = ws0*sb_tests0;
      ws1 = ws1*sb_tests1;
      ws2 = ws2*sb_tests2;
      ws3 = ws3*sb_tests3;
      ws4 = ws4*sb_tests4;
      ws5 = ws5*sb_tests5;
      ws6 = ws6*sb_tests6;
      ws7 = ws7*sb_tests7;
      ws8 = ws8*sb_tests8;
      ws9 = ws9*sb_tests9;
      wsa = wsa*sb_testsa;
      wsb = wsb*sb_testsb;
      wsc = wsc*sb_testsc;
      wsd = wsd*sb_testsd;
      wse = wse*sb_testse;
      wsf = wsf*sb_testsf;

      U_1[k] = ws0+ws1+ws2+ws3
	       +ws4+ws5+ws6+ws7
	       +ws8+ws9+wsa+wsb
	       +wsc+wsd+wse+wsf;

      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_stag-1 & N-1);
      ind_2y = (i_2_stag & N-1);
      ind_2z = (i_2_stag+1 & N-1);
      ind_2w = (i_2_stag+2 & N-1);

      sb_tests0 = u_2[ind_1x+N*ind_2x];
      sb_tests1 = u_2[ind_1y+N*ind_2x];
      sb_tests2 = u_2[ind_1z+N*ind_2x];
      sb_tests3 = u_2[ind_1w+N*ind_2x];
      sb_tests4 = u_2[ind_1x+N*ind_2y];
      sb_tests5 = u_2[ind_1y+N*ind_2y];
      sb_tests6 = u_2[ind_1z+N*ind_2y];
      sb_tests7 = u_2[ind_1w+N*ind_2y];
      sb_tests8 = u_2[ind_1x+N*ind_2z];
      sb_tests9 = u_2[ind_1y+N*ind_2z];
      sb_testsa = u_2[ind_1z+N*ind_2z];
      sb_testsb = u_2[ind_1w+N*ind_2z];
      sb_testsc = u_2[ind_1x+N*ind_2w];
      sb_testsd = u_2[ind_1y+N*ind_2w];
      sb_testse = u_2[ind_1z+N*ind_2w];
      sb_testsf = u_2[ind_1w+N*ind_2w];

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      ws0 = ws0*sb_tests0;
      ws1 = ws1*sb_tests1;
      ws2 = ws2*sb_tests2;
      ws3 = ws3*sb_tests3;
      ws4 = ws4*sb_tests4;
      ws5 = ws5*sb_tests5;
      ws6 = ws6*sb_tests6;
      ws7 = ws7*sb_tests7;
      ws8 = ws8*sb_tests8;
      ws9 = ws9*sb_tests9;
      wsa = wsa*sb_testsa;
      wsb = wsb*sb_testsb;
      wsc = wsc*sb_testsc;
      wsd = wsd*sb_testsd;
      wse = wse*sb_testse;
      wsf = wsf*sb_testsf;

      U_2[k] = ws0+ws1+ws2+ws3
	       +ws4+ws5+ws6+ws7
	       +ws8+ws9+wsa+wsb
	       +wsc+wsd+wse+wsf;
    }
    __kernel void spread(__global const int *firstn1, __global const int *nextn1, __global const int *firstn2, __global const int *nextn2, __global const double *F_1, __global const double *F_2, __global const double *Mu, __global const double *Rho, __global const double *X_1, __global const double *X_2, const unsigned N, const double dx, const double dtheta, const unsigned shift0, const unsigned shift1, __global double *f_1, __global double *f_2, __global double *mu_1, __global double *mu_2, __global double *rho_1, __global double *rho_2)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      double c = dtheta/(dx*dx);
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;

      double f_10, f_11, f_12, f_13;
      double f_14, f_15, f_16, f_17;
      double f_18, f_19, f_1a, f_1b;
      double f_1c, f_1d, f_1e, f_1f;

      double f_20, f_21, f_22, f_23;
      double f_24, f_25, f_26, f_27;
      double f_28, f_29, f_2a, f_2b;
      double f_2c, f_2d, f_2e, f_2f;

      double mu_10, mu_11, mu_12, mu_13;
      double mu_14, mu_15, mu_16, mu_17;
      double mu_18, mu_19, mu_1a, mu_1b;
      double mu_1c, mu_1d, mu_1e, mu_1f;

      double mu_20, mu_21, mu_22, mu_23;
      double mu_24, mu_25, mu_26, mu_27;
      double mu_28, mu_29, mu_2a, mu_2b;
      double mu_2c, mu_2d, mu_2e, mu_2f;

      double rho_10, rho_11, rho_12, rho_13;
      double rho_14, rho_15, rho_16, rho_17;
      double rho_18, rho_19, rho_1a, rho_1b;
      double rho_1c, rho_1d, rho_1e, rho_1f;

      double rho_20, rho_21, rho_22, rho_23;
      double rho_24, rho_25, rho_26, rho_27;
      double rho_28, rho_29, rho_2a, rho_2b;
      double rho_2c, rho_2d, rho_2e, rho_2f;

      double v1x, v1y, v1z, v1w;
      double v2x, v2y, v2z, v2w;

      double q1, q2;

      double ws0, ws1, ws2, ws3;
      double ws4, ws5, ws6, ws7;
      double ws8, ws9, wsa, wsb;
      double wsc, wsd, wse, wsf;

      if (k1 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_10=f_1[ind_1x+N*ind_2x];
      f_11=f_1[ind_1y+N*ind_2x];
      f_12=f_1[ind_1z+N*ind_2x];
      f_13=f_1[ind_1w+N*ind_2x];
      f_14=f_1[ind_1x+N*ind_2y];
      f_15=f_1[ind_1y+N*ind_2y];
      f_16=f_1[ind_1z+N*ind_2y];
      f_17=f_1[ind_1w+N*ind_2y];
      f_18=f_1[ind_1x+N*ind_2z];
      f_19=f_1[ind_1y+N*ind_2z];
      f_1a=f_1[ind_1z+N*ind_2z];
      f_1b=f_1[ind_1w+N*ind_2z];
      f_1c=f_1[ind_1x+N*ind_2w];
      f_1d=f_1[ind_1y+N*ind_2w];
      f_1e=f_1[ind_1z+N*ind_2w];
      f_1f=f_1[ind_1w+N*ind_2w];

      mu_10=mu_1[ind_1x+N*ind_2x];
      mu_11=mu_1[ind_1y+N*ind_2x];
      mu_12=mu_1[ind_1z+N*ind_2x];
      mu_13=mu_1[ind_1w+N*ind_2x];
      mu_14=mu_1[ind_1x+N*ind_2y];
      mu_15=mu_1[ind_1y+N*ind_2y];
      mu_16=mu_1[ind_1z+N*ind_2y];
      mu_17=mu_1[ind_1w+N*ind_2y];
      mu_18=mu_1[ind_1x+N*ind_2z];
      mu_19=mu_1[ind_1y+N*ind_2z];
      mu_1a=mu_1[ind_1z+N*ind_2z];
      mu_1b=mu_1[ind_1w+N*ind_2z];
      mu_1c=mu_1[ind_1x+N*ind_2w];
      mu_1d=mu_1[ind_1y+N*ind_2w];
      mu_1e=mu_1[ind_1z+N*ind_2w];
      mu_1f=mu_1[ind_1w+N*ind_2w];

      rho_10=rho_1[ind_1x+N*ind_2x];
      rho_11=rho_1[ind_1y+N*ind_2x];
      rho_12=rho_1[ind_1z+N*ind_2x];
      rho_13=rho_1[ind_1w+N*ind_2x];
      rho_14=rho_1[ind_1x+N*ind_2y];
      rho_15=rho_1[ind_1y+N*ind_2y];
      rho_16=rho_1[ind_1z+N*ind_2y];
      rho_17=rho_1[ind_1w+N*ind_2y];
      rho_18=rho_1[ind_1x+N*ind_2z];
      rho_19=rho_1[ind_1y+N*ind_2z];
      rho_1a=rho_1[ind_1z+N*ind_2z];
      rho_1b=rho_1[ind_1w+N*ind_2z];
      rho_1c=rho_1[ind_1x+N*ind_2w];
      rho_1d=rho_1[ind_1y+N*ind_2w];
      rho_1e=rho_1[ind_1z+N*ind_2w];
      rho_1f=rho_1[ind_1w+N*ind_2w];
      }

      while(k1 != -1)
      {
      double s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      double r_1_stag = s_1_stag-i_1_stag;
      double s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
      double r_2_c = s_2_c-i_2_c;
      
      q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      double my_F_1 = F_1[k1];
      f_10 += c*my_F_1*ws0;
      f_11 += c*my_F_1*ws1;
      f_12 += c*my_F_1*ws2;
      f_13 += c*my_F_1*ws3;
      f_14 += c*my_F_1*ws4;
      f_15 += c*my_F_1*ws5;
      f_16 += c*my_F_1*ws6;
      f_17 += c*my_F_1*ws7;
      f_18 += c*my_F_1*ws8;
      f_19 += c*my_F_1*ws9;
      f_1a += c*my_F_1*wsa;
      f_1b += c*my_F_1*wsb;
      f_1c += c*my_F_1*wsc;
      f_1d += c*my_F_1*wsd;
      f_1e += c*my_F_1*wse;
      f_1f += c*my_F_1*wsf;

      double my_Mu = Mu[k1];
      mu_10 += c*my_Mu*ws0;
      mu_11 += c*my_Mu*ws1;
      mu_12 += c*my_Mu*ws2;
      mu_13 += c*my_Mu*ws3;
      mu_14 += c*my_Mu*ws4;
      mu_15 += c*my_Mu*ws5;
      mu_16 += c*my_Mu*ws6;
      mu_17 += c*my_Mu*ws7;
      mu_18 += c*my_Mu*ws8;
      mu_19 += c*my_Mu*ws9;
      mu_1a += c*my_Mu*wsa;
      mu_1b += c*my_Mu*wsb;
      mu_1c += c*my_Mu*wsc;
      mu_1d += c*my_Mu*wsd;
      mu_1e += c*my_Mu*wse;
      mu_1f += c*my_Mu*wsf;

      double my_Rho = Rho[k1];
      rho_10 += c*my_Rho*ws0;
      rho_11 += c*my_Rho*ws1;
      rho_12 += c*my_Rho*ws2;
      rho_13 += c*my_Rho*ws3;
      rho_14 += c*my_Rho*ws4;
      rho_15 += c*my_Rho*ws5;
      rho_16 += c*my_Rho*ws6;
      rho_17 += c*my_Rho*ws7;
      rho_18 += c*my_Rho*ws8;
      rho_19 += c*my_Rho*ws9;
      rho_1a += c*my_Rho*wsa;
      rho_1b += c*my_Rho*wsb;
      rho_1c += c*my_Rho*wsc;
      rho_1d += c*my_Rho*wsd;
      rho_1e += c*my_Rho*wse;
      rho_1f += c*my_Rho*wsf;

      k1 = nextn1[k1];
      if (k1 == -1){
      f_1[ind_1x+N*ind_2x]=f_10;
      f_1[ind_1y+N*ind_2x]=f_11;
      f_1[ind_1z+N*ind_2x]=f_12;
      f_1[ind_1w+N*ind_2x]=f_13;
      f_1[ind_1x+N*ind_2y]=f_14;
      f_1[ind_1y+N*ind_2y]=f_15;
      f_1[ind_1z+N*ind_2y]=f_16;
      f_1[ind_1w+N*ind_2y]=f_17;
      f_1[ind_1x+N*ind_2z]=f_18;
      f_1[ind_1y+N*ind_2z]=f_19;
      f_1[ind_1z+N*ind_2z]=f_1a;
      f_1[ind_1w+N*ind_2z]=f_1b;
      f_1[ind_1x+N*ind_2w]=f_1c;
      f_1[ind_1y+N*ind_2w]=f_1d;
      f_1[ind_1z+N*ind_2w]=f_1e;
      f_1[ind_1w+N*ind_2w]=f_1f;

      mu_1[ind_1x+N*ind_2x]=mu_10;
      mu_1[ind_1y+N*ind_2x]=mu_11;
      mu_1[ind_1z+N*ind_2x]=mu_12;
      mu_1[ind_1w+N*ind_2x]=mu_13;
      mu_1[ind_1x+N*ind_2y]=mu_14;
      mu_1[ind_1y+N*ind_2y]=mu_15;
      mu_1[ind_1z+N*ind_2y]=mu_16;
      mu_1[ind_1w+N*ind_2y]=mu_17;
      mu_1[ind_1x+N*ind_2z]=mu_18;
      mu_1[ind_1y+N*ind_2z]=mu_19;
      mu_1[ind_1z+N*ind_2z]=mu_1a;
      mu_1[ind_1w+N*ind_2z]=mu_1b;
      mu_1[ind_1x+N*ind_2w]=mu_1c;
      mu_1[ind_1y+N*ind_2w]=mu_1d;
      mu_1[ind_1z+N*ind_2w]=mu_1e;
      mu_1[ind_1w+N*ind_2w]=mu_1f;

      rho_1[ind_1x+N*ind_2x]=rho_10;
      rho_1[ind_1y+N*ind_2x]=rho_11;
      rho_1[ind_1z+N*ind_2x]=rho_12;
      rho_1[ind_1w+N*ind_2x]=rho_13;
      rho_1[ind_1x+N*ind_2y]=rho_14;
      rho_1[ind_1y+N*ind_2y]=rho_15;
      rho_1[ind_1z+N*ind_2y]=rho_16;
      rho_1[ind_1w+N*ind_2y]=rho_17;
      rho_1[ind_1x+N*ind_2z]=rho_18;
      rho_1[ind_1y+N*ind_2z]=rho_19;
      rho_1[ind_1z+N*ind_2z]=rho_1a;
      rho_1[ind_1w+N*ind_2z]=rho_1b;
      rho_1[ind_1x+N*ind_2w]=rho_1c;
      rho_1[ind_1y+N*ind_2w]=rho_1d;
      rho_1[ind_1z+N*ind_2w]=rho_1e;
      rho_1[ind_1w+N*ind_2w]=rho_1f;
      }
      }
      
      if (k2 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_20=f_2[ind_1x+N*ind_2x];
      f_21=f_2[ind_1y+N*ind_2x];
      f_22=f_2[ind_1z+N*ind_2x];
      f_23=f_2[ind_1w+N*ind_2x];
      f_24=f_2[ind_1x+N*ind_2y];
      f_25=f_2[ind_1y+N*ind_2y];
      f_26=f_2[ind_1z+N*ind_2y];
      f_27=f_2[ind_1w+N*ind_2y];
      f_28=f_2[ind_1x+N*ind_2z];
      f_29=f_2[ind_1y+N*ind_2z];
      f_2a=f_2[ind_1z+N*ind_2z];
      f_2b=f_2[ind_1w+N*ind_2z];
      f_2c=f_2[ind_1x+N*ind_2w];
      f_2d=f_2[ind_1y+N*ind_2w];
      f_2e=f_2[ind_1z+N*ind_2w];
      f_2f=f_2[ind_1w+N*ind_2w];

      mu_20=mu_2[ind_1x+N*ind_2x];
      mu_21=mu_2[ind_1y+N*ind_2x];
      mu_22=mu_2[ind_1z+N*ind_2x];
      mu_23=mu_2[ind_1w+N*ind_2x];
      mu_24=mu_2[ind_1x+N*ind_2y];
      mu_25=mu_2[ind_1y+N*ind_2y];
      mu_26=mu_2[ind_1z+N*ind_2y];
      mu_27=mu_2[ind_1w+N*ind_2y];
      mu_28=mu_2[ind_1x+N*ind_2z];
      mu_29=mu_2[ind_1y+N*ind_2z];
      mu_2a=mu_2[ind_1z+N*ind_2z];
      mu_2b=mu_2[ind_1w+N*ind_2z];
      mu_2c=mu_2[ind_1x+N*ind_2w];
      mu_2d=mu_2[ind_1y+N*ind_2w];
      mu_2e=mu_2[ind_1z+N*ind_2w];
      mu_2f=mu_2[ind_1w+N*ind_2w];

      rho_20=rho_2[ind_1x+N*ind_2x];
      rho_21=rho_2[ind_1y+N*ind_2x];
      rho_22=rho_2[ind_1z+N*ind_2x];
      rho_23=rho_2[ind_1w+N*ind_2x];
      rho_24=rho_2[ind_1x+N*ind_2y];
      rho_25=rho_2[ind_1y+N*ind_2y];
      rho_26=rho_2[ind_1z+N*ind_2y];
      rho_27=rho_2[ind_1w+N*ind_2y];
      rho_28=rho_2[ind_1x+N*ind_2z];
      rho_29=rho_2[ind_1y+N*ind_2z];
      rho_2a=rho_2[ind_1z+N*ind_2z];
      rho_2b=rho_2[ind_1w+N*ind_2z];
      rho_2c=rho_2[ind_1x+N*ind_2w];
      rho_2d=rho_2[ind_1y+N*ind_2w];
      rho_2e=rho_2[ind_1z+N*ind_2w];
      rho_2f=rho_2[ind_1w+N*ind_2w];
      }
      
      while (k2 != -1){
      double s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      double r_1_c = s_1_c-i_1_c;
      double s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      double r_2_stag = s_2_stag-i_2_stag;

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      double my_F_2 = F_2[k2];
      f_20 += c*my_F_2*ws0;
      f_21 += c*my_F_2*ws1;
      f_22 += c*my_F_2*ws2;
      f_23 += c*my_F_2*ws3;
      f_24 += c*my_F_2*ws4;
      f_25 += c*my_F_2*ws5;
      f_26 += c*my_F_2*ws6;
      f_27 += c*my_F_2*ws7;
      f_28 += c*my_F_2*ws8;
      f_29 += c*my_F_2*ws9;
      f_2a += c*my_F_2*wsa;
      f_2b += c*my_F_2*wsb;
      f_2c += c*my_F_2*wsc;
      f_2d += c*my_F_2*wsd;
      f_2e += c*my_F_2*wse;
      f_2f += c*my_F_2*wsf;

      double my_Mu = Mu[k2];
      mu_20 += c*my_Mu*ws0;
      mu_21 += c*my_Mu*ws1;
      mu_22 += c*my_Mu*ws2;
      mu_23 += c*my_Mu*ws3;
      mu_24 += c*my_Mu*ws4;
      mu_25 += c*my_Mu*ws5;
      mu_26 += c*my_Mu*ws6;
      mu_27 += c*my_Mu*ws7;
      mu_28 += c*my_Mu*ws8;
      mu_29 += c*my_Mu*ws9;
      mu_2a += c*my_Mu*wsa;
      mu_2b += c*my_Mu*wsb;
      mu_2c += c*my_Mu*wsc;
      mu_2d += c*my_Mu*wsd;
      mu_2e += c*my_Mu*wse;
      mu_2f += c*my_Mu*wsf;

      double my_Rho = Rho[k2];
      rho_20 += c*my_Rho*ws0;
      rho_21 += c*my_Rho*ws1;
      rho_22 += c*my_Rho*ws2;
      rho_23 += c*my_Rho*ws3;
      rho_24 += c*my_Rho*ws4;
      rho_25 += c*my_Rho*ws5;
      rho_26 += c*my_Rho*ws6;
      rho_27 += c*my_Rho*ws7;
      rho_28 += c*my_Rho*ws8;
      rho_29 += c*my_Rho*ws9;
      rho_2a += c*my_Rho*wsa;
      rho_2b += c*my_Rho*wsb;
      rho_2c += c*my_Rho*wsc;
      rho_2d += c*my_Rho*wsd;
      rho_2e += c*my_Rho*wse;
      rho_2f += c*my_Rho*wsf;

      k2 = nextn2[k2];
      if (k2 == -1){
      f_2[ind_1x+N*ind_2x]=f_20;
      f_2[ind_1y+N*ind_2x]=f_21;
      f_2[ind_1z+N*ind_2x]=f_22;
      f_2[ind_1w+N*ind_2x]=f_23;
      f_2[ind_1x+N*ind_2y]=f_24;
      f_2[ind_1y+N*ind_2y]=f_25;
      f_2[ind_1z+N*ind_2y]=f_26;
      f_2[ind_1w+N*ind_2y]=f_27;
      f_2[ind_1x+N*ind_2z]=f_28;
      f_2[ind_1y+N*ind_2z]=f_29;
      f_2[ind_1z+N*ind_2z]=f_2a;
      f_2[ind_1w+N*ind_2z]=f_2b;
      f_2[ind_1x+N*ind_2w]=f_2c;
      f_2[ind_1y+N*ind_2w]=f_2d;
      f_2[ind_1z+N*ind_2w]=f_2e;
      f_2[ind_1w+N*ind_2w]=f_2f;

      mu_2[ind_1x+N*ind_2x]=mu_20;
      mu_2[ind_1y+N*ind_2x]=mu_21;
      mu_2[ind_1z+N*ind_2x]=mu_22;
      mu_2[ind_1w+N*ind_2x]=mu_23;
      mu_2[ind_1x+N*ind_2y]=mu_24;
      mu_2[ind_1y+N*ind_2y]=mu_25;
      mu_2[ind_1z+N*ind_2y]=mu_26;
      mu_2[ind_1w+N*ind_2y]=mu_27;
      mu_2[ind_1x+N*ind_2z]=mu_28;
      mu_2[ind_1y+N*ind_2z]=mu_29;
      mu_2[ind_1z+N*ind_2z]=mu_2a;
      mu_2[ind_1w+N*ind_2z]=mu_2b;
      mu_2[ind_1x+N*ind_2w]=mu_2c;
      mu_2[ind_1y+N*ind_2w]=mu_2d;
      mu_2[ind_1z+N*ind_2w]=mu_2e;
      mu_2[ind_1w+N*ind_2w]=mu_2f;

      rho_2[ind_1x+N*ind_2x]=rho_20;
      rho_2[ind_1y+N*ind_2x]=rho_21;
      rho_2[ind_1z+N*ind_2x]=rho_22;
      rho_2[ind_1w+N*ind_2x]=rho_23;
      rho_2[ind_1x+N*ind_2y]=rho_24;
      rho_2[ind_1y+N*ind_2y]=rho_25;
      rho_2[ind_1z+N*ind_2y]=rho_26;
      rho_2[ind_1w+N*ind_2y]=rho_27;
      rho_2[ind_1x+N*ind_2z]=rho_28;
      rho_2[ind_1y+N*ind_2z]=rho_29;
      rho_2[ind_1z+N*ind_2z]=rho_2a;
      rho_2[ind_1w+N*ind_2z]=rho_2b;
      rho_2[ind_1x+N*ind_2w]=rho_2c;
      rho_2[ind_1y+N*ind_2w]=rho_2d;
      rho_2[ind_1z+N*ind_2w]=rho_2e;
      rho_2[ind_1w+N*ind_2w]=rho_2f;
      }
      }
    }
    __kernel void updaten(__global int *firstn1, __global int *nextn1, __global int *firstn2, __global int *nextn2, __global const double *X_1, __global const double *X_2, const unsigned N, const double dx, const unsigned shift0, const unsigned shift1)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      uint ind_1, ind_2;
      int nextnold;

      int nprev1 = -1;
      while(k1 != -1)
      {
      double s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      double s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
       
      ind_1 = (i_1_stag & N-1);
      ind_2 = (i_2_c & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev1 = k1;
      k1 = nextn1[k1];
      }
      else {
      nextnold = nextn1[k1];
      nextn1[k1] = firstn1[ind_1+N*ind_2];
      firstn1[ind_1+N*ind_2] = k1;
      k1 = nextnold;
          if (nprev1 == -1) {
          firstn1[i+N*j] = nextnold;
          }
          else {
          nextn1[nprev1] = nextnold;
          }
      }
      }

      int nprev2 = -1;
      while(k2 != -1)
      {
      double s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      double s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      
      ind_1 = (i_1_c & N-1);
      ind_2 = (i_2_stag & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev2 = k2;
      k2 = nextn2[k2];
      }
      else {
      nextnold = nextn2[k2];
      nextn2[k2] = firstn2[ind_1+N*ind_2];
      firstn2[ind_1+N*ind_2] = k2;
      k2 = nextnold;
          if (nprev2 == -1) {
          firstn2[i+N*j] = nextnold;
          }
          else {
          nextn2[nprev2] = nextnold;
          }
      }
      }
    }
    __kernel void calc_lagF(__global const double *X_1, __global const double *X_2, __global const double *X_1_teth, __global const double *X_2_teth, __global const unsigned *i_teth, __global const unsigned *i_t, __global const unsigned *ip_t, __global const unsigned *im_t, const double k_teth, const double k_s, const double k_b, const double ds, __global double *F_1, __global double *F_2) 
    {
    unsigned i = get_global_id(0);
    unsigned N_r = get_global_size(0);
    N_r = N_r/2;

    F_1[i] = -k_teth/ds*(X_1[i_teth[i]]-X_1_teth[i]) 
             -k_s*((sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)         
             *(X_1[i_t[i]]-X_1[im_t[i]])/sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))
             -(sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_1[ip_t[i]]-X_1[i_t[i]])/sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]])))
             -k_b*(X_1[ip_t[ip_t[i]]]-4*X_1[ip_t[i]]+X_1[im_t[im_t[i]]]-4*X_1[im_t[i]]+6*X_1[i_t[i]])/(ds*ds*ds);
    F_2[i] = -k_teth/ds*(X_2[i_teth[i]]-X_2_teth[i]) 
             -k_s*((sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)         
             *(X_2[i_t[i]]-X_2[im_t[i]])/sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))
             -(sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_2[ip_t[i]]-X_2[i_t[i]])/sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]])))
             -k_b*(X_2[ip_t[ip_t[i]]]-4*X_2[ip_t[i]]+X_2[im_t[im_t[i]]]-4*X_2[im_t[i]]+6*X_2[i_t[i]])/(ds*ds*ds);

    if (i_t[i] == 1 || i_t[i] == 1+N_r){
    F_1[i] += k_b/(ds*ds*ds);
    }
    else if (i_t[i] == N_r-2 || i_t[i] == 2*N_r-2){
    F_1[i] += -k_b/(ds*ds*ds);
    }
    else if (i_t[i] == 0 || i_t[i] == N_r){
    F_1[i] += -3.0*k_b/(ds*ds*ds)
             +k_s*(sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)
             *(X_1[i_t[i]]-X_1[im_t[i]])/sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))
             -k_s*(sqrt((X_1[i_t[i]]-X_1[im_t[i]]+1.0)*(X_1[i_t[i]]-X_1[im_t[i]]+1.0)+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)
             *(X_1[i_t[i]]-X_1[im_t[i]]+1.0)/sqrt((X_1[i_t[i]]-X_1[im_t[i]]+1.0)*(X_1[i_t[i]]-X_1[im_t[i]]+1.0)+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]));
    F_2[i] += k_s*(sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)
             *(X_2[i_t[i]]-X_2[im_t[i]])/sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))
              -k_s*(sqrt((X_1[i_t[i]]-X_1[im_t[i]]+1.0)*(X_1[i_t[i]]-X_1[im_t[i]]+1.0)+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)
             *(X_2[i_t[i]]-X_2[im_t[i]])/sqrt((X_1[i_t[i]]-X_1[im_t[i]]+1.0)*(X_1[i_t[i]]-X_1[im_t[i]]+1.0)+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]));
    }
    else if (i_t[i] == N_r-1 || i_t[i] == 2*N_r-1){
    F_1[i] += 3.0*k_b/(ds*ds*ds)
             -k_s*(sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_1[ip_t[i]]-X_1[i_t[i]])/sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))
             +k_s*(sqrt((X_1[ip_t[i]]+1.0-X_1[i_t[i]])*(X_1[ip_t[i]]+1.0-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_1[ip_t[i]]+1.0-X_1[i_t[i]])/sqrt((X_1[ip_t[i]]+1.0-X_1[i_t[i]])*(X_1[ip_t[i]]+1.0-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]));
    F_2[i] += -k_s*(sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_2[ip_t[i]]-X_2[i_t[i]])/sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))
             +k_s*(sqrt((X_1[ip_t[i]]+1.0-X_1[i_t[i]])*(X_1[ip_t[i]]+1.0-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_2[ip_t[i]]-X_2[i_t[i]])/sqrt((X_1[ip_t[i]]+1.0-X_1[i_t[i]])*(X_1[ip_t[i]]+1.0-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]));
    }
    }
    """).build()

solve_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    #pragma OPENCL EXTENSION cl_khr_fp64 : enable
    //#pragma OPENCL EXTENSION cl_amd_fp64 : enable
    #define complex_ctr(x, y) (double2)(x, y)
    #define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
    #define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
    __kernel void solve(__global const double2 *u_1_tf, __global const double2 *u_2_tf, __global const double2 *L_tf,
    __global const double2 *f_1_tf, __global const double2 *f_2_tf,
    const double dt, const double rho, const double mu, const unsigned N,
    __global double2 *u_1_tf_tilde, __global double2 *u_2_tf_tilde)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned base = i + N*j;
      double2 m;
      m.x = 1.0+mu*dt/(2.0*rho)*L_tf[base].x; 
      m.y = mu*dt/(2.0*rho)*L_tf[base].y; 
      double2 rhs_1 = complex_mul(m, u_1_tf[base])+dt/rho*f_1_tf[base];
      double2 rhs_2 = complex_mul(m, u_2_tf[base])+dt/rho*f_2_tf[base];
      double2 lhs;
      lhs.x = 1.0-mu*dt/(2.0*rho)*L_tf[base].x;
      lhs.y = -mu*dt/(2.0*rho)*L_tf[base].y;
      u_1_tf_tilde[base] = complex_div(rhs_1, lhs);
      u_2_tf_tilde[base] = complex_div(rhs_2, lhs);
    }
    __kernel void dal_term(__global const double *u_1, __global const double *u_2, __global const double *u_1_o, __global const double *u_2_o, __global const double *ss_1, __global const double *ss_2, __global const double *rho_1, __global const double *rho_2,
    const double dt, const unsigned N,
    __global double *f_1, __global double *f_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned base = i + N*j;
      f_1[base] = rho_1[base]*((u_1[base]-u_1_o[base])/dt+ss_1[base]);
      f_2[base] = rho_2[base]*((u_2[base]-u_2_o[base])/dt+ss_2[base]);
    }
    __kernel void vc_laplace(__global const double *u_1,
    __global const double *u_2,
    __global const double *mu_1,
    __global const double *mu_2, const double dx, const unsigned N,
    __global double *V_1, __global double *V_2)
    {
    unsigned i = get_global_id(0);
    unsigned j = get_global_id(1);
    unsigned base = i + N*j;
    unsigned base1 = (i+1 & N-1) + N*j;
    unsigned base2 = (i-1 & N-1) + N*j;
    unsigned base3 = i + N*(j+1 & N-1);
    unsigned base4 = i + N*(j-1 & N-1);
    unsigned base5 = (i-1 & N-1) + N*(j+1 & N-1);
    unsigned base6 = (i+1 & N-1) + N*(j-1 & N-1);

    V_1[base] = 1.0/(dx*dx)*(2.0*(0.5*(mu_1[base]+mu_1[base1])*(u_1[base1]-u_1[base])
               -0.5*(mu_1[base]+mu_1[base2])*(u_1[base]-u_1[base2]))
               +0.5*(mu_1[base3]+mu_1[base])*(u_1[base3]-u_1[base])
               -0.5*(mu_1[base]+mu_1[base4])*(u_1[base]-u_1[base4]) 
               +0.5*(mu_1[base3]+mu_1[base])*(u_2[base3]-u_2[base5])
               -0.5*(mu_1[base]+mu_1[base4])*(u_2[base]-u_2[base2]));

    V_2[base] = 1.0/(dx*dx)*(0.5*(mu_2[base]+mu_2[base1])*(u_2[base1]-u_2[base])
               -0.5*(mu_2[base]+mu_2[base2])*(u_2[base]-u_2[base2])
               +0.5*(mu_2[base]+mu_2[base1])*(u_1[base1]-u_1[base6])
               -0.5*(mu_2[base]+mu_2[base2])*(u_1[base]-u_1[base4])
               +2*(0.5*(mu_2[base]+mu_2[base3])*(u_2[base3]-u_2[base])
               -0.5*(mu_2[base]+mu_2[base4])*(u_2[base]-u_2[base4])));

    }   
    """).build()

prg = cl.Program(ctx, """
    #pragma OPENCL EXTENSION cl_khr_fp64 : enable
    //#pragma OPENCL EXTENSION cl_amd_fp64 : enable
    __kernel void p_filt(__global double2 *u)
    {
      u[0].x = 0;
      u[0].y = 0;
    }
    __kernel void L_adj_1(__global double2 *L_tf)
    {
      L_tf[0].x = 0.0;
      L_tf[0].y = 0.0;
    }
    __kernel void L_adj_2(__global double2 *L_tf)
    {
      L_tf[0].x = 1.0;
      L_tf[0].y = 0.0;
    }
    __kernel void nyq_filt(__global double2 *u, const unsigned N)
    {
      unsigned i = get_global_id(0);
      unsigned base = N/2 + N*i; 
      unsigned base1 = i + N*N/2; 
   
      u[base].x = 0;
      u[base].y = 0;
      u[base1].x = 0;
      u[base1].y = 0;
    }
    __kernel void vort(__global const double *u_1,
    __global const double *u_2,
    const double dx, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned base = i + N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = i+N*(j-1 & N-1);
      f[base] = 1/dx*(u_2[base]-u_2[base1]-(u_1[base]-u_1[base2]));
    }
    __kernel void div(__global const double *u_1,
    __global const double *u_2,
    const double dx, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //f[i,j] = 1/dx*(u_1[i_2[i,j],ip_1[i,j]]-u_1[i,j]+u_2[ip_2[i,j],i_1[i,j]]-u_2[i,j];
      unsigned base = i + N*j;
      unsigned base1 = (i+1 & N-1)+N*j;
      unsigned base2 = i+N*(j+1 & N-1);
      f[base] = 1/dx*(u_1[base1]-u_1[base]+u_2[base2]-u_2[base]);
    }
    __kernel void grad(__global const double *phi,
    const double dx, const unsigned N,
    __global double *g_1, __global double *g_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //grad_1 = 1/dx*(phi[i_2,i_1]-phi[i_2,im_1])
      //grad_2 = 1/dx*(phi[i_2,i_1]-phi[im_2,i_1])
      unsigned base = i+N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = i+N*(j-1 & N-1);
      g_1[base] = 1/dx*(phi[base]-phi[base1]);
      g_2[base] = 1/dx*(phi[base]-phi[base2]);
    }
    __kernel void ss_conv(__global const double *u_1,
    __global const double *u_2, const double dx, const unsigned N,
    __global double *S_1, __global double *S_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //Compute S = 0.5*(div(UU) + (U*grad)U).
      unsigned base = i+N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = (i+1 & N-1)+N*j;
      unsigned base3 = i+N*(j-1 & N-1);
      unsigned base4 = i+N*(j+1 & N-1);
      unsigned base5 = (i-1 & N-1)+N*(j+1 & N-1);
      unsigned base6 = (i+1 & N-1)+N*(j-1 & N-1);
      double D_x, D_y, A_x, A_y, D, A;
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]+u_1[base]) -
          (u_1[base]+u_1[base1])*(u_1[base]+u_1[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base5])*(u_1[base4]+u_1[base]) -
          (u_2[base]+u_2[base1])*(u_1[base]+u_1[base3]));
      D = D_x + D_y;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]-u_1[base]) +
          (u_1[base]+u_1[base1])*(u_1[base]-u_1[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base5])*(u_1[base4]-u_1[base]) +
            (u_2[base]+u_2[base1])*(u_1[base]-u_1[base3]));
      A = A_x + A_y;
      S_1[base] = 0.5*(D+A);
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]+u_2[base]) -
          (u_1[base]+u_1[base3])*(u_2[base]+u_2[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base])*(u_2[base4]+u_2[base]) -
          (u_2[base]+u_2[base3])*(u_2[base]+u_2[base3]));
      D = D_x + D_y;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]-u_2[base]) +
          (u_1[base]+u_1[base3])*(u_2[base]-u_2[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base])*(u_2[base4]-u_2[base]) +
            (u_2[base]+u_2[base3])*(u_2[base]-u_2[base3]));
      A = A_x + A_y;
      S_2[base] = 0.5*(D+A);
    }
    __kernel void interpctolr(__global const double *phi, __global const unsigned *i_2,  __global const unsigned *im_1, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[i_2,im_1])

      unsigned base = i+N*j;
      unsigned base1 = im_1[base]+N*i_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interpctotb(__global const double *phi, __global const unsigned *i_1, __global const unsigned *im_2, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[im_2,i_1])

      unsigned base = i+N*j;
      unsigned base1 = i_1[base]+N*im_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interplrtoc(__global const double *phi, __global const unsigned *i_2, __global const unsigned *ip_1, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[i_2,ip_1])

      unsigned base = i+N*j;
      unsigned base1 = ip_1[base]+N*i_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interptbtoc(__global const double *phi, __global const unsigned *i_1,__global const unsigned *ip_2, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[ip_2,i_1])

      unsigned base = i+N*j;
      unsigned base1 = i_1[base]+N*ip_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    """).build()

calc_lagF_kernel = ib_prg.calc_lagF
calc_lagF_kernel.set_scalar_arg_dtypes([None]*8 + [np.float64]*4 + [None]*2)

div_kernel = prg.div
div_kernel.set_scalar_arg_dtypes([None]*2 + [np.float64, np.uint32] + [None])

vort_kernel = prg.vort
vort_kernel.set_scalar_arg_dtypes([None]*2 + [np.float64, np.uint32] + [None])

L_adj_1_kernel = prg.L_adj_1
L_adj_1_kernel.set_scalar_arg_dtypes([None])

L_adj_2_kernel = prg.L_adj_2
L_adj_2_kernel.set_scalar_arg_dtypes([None])

p_filt_kernel = prg.p_filt
p_filt_kernel.set_scalar_arg_dtypes([None])

nyq_filt_kernel = prg.nyq_filt
nyq_filt_kernel.set_scalar_arg_dtypes([None] + [np.uint32])

solve_kernel = solve_prg.solve
solve_kernel.set_scalar_arg_dtypes([None]*5 + [np.float64]*3 + [np.uint32] + [None]*2)

vc_laplace_kernel = solve_prg.vc_laplace
vc_laplace_kernel.set_scalar_arg_dtypes([None]*4 + [np.float64] + [np.uint32] + [None]*2)

dal_term_kernel = solve_prg.dal_term
dal_term_kernel.set_scalar_arg_dtypes([None]*8 + [np.float64] + [np.uint32] + [None]*2)

grad_kernel = prg.grad
grad_kernel.set_scalar_arg_dtypes([None]*1 + [np.float64, np.uint32] + [None]*2)

ss_conv_kernel = prg.ss_conv
ss_conv_kernel.set_scalar_arg_dtypes([None]*2 + [np.float64, np.uint32] + [None]*2)

interp_kernel = ib_prg.interp
interp_kernel.set_scalar_arg_dtypes([None]*4 + [np.uint32, np.float64] + [None]*2)
spread_kernel = ib_prg.spread
spread_kernel.set_scalar_arg_dtypes([None]*10 + [np.uint32] + [np.float64]*2 + [np.uint32]*2 + [None]*6)
updaten_kernel = ib_prg.updaten
updaten_kernel.set_scalar_arg_dtypes([None]*6 + [np.uint32] + [np.float64] + [np.uint32]*2)

interpctolr_kernel = prg.interpctolr
interpctolr_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interpctotb_kernel = prg.interpctotb
interpctotb_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interplrtoc_kernel = prg.interplrtoc
interplrtoc_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interptbtoc_kernel = prg.interptbtoc
interptbtoc_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])

for ii in range(start_size,end_size+1): 
    #set up grid
    N = 2**ii
    dx = 1.0/N
    rho_0 = 1.0

    x_1_stag = np.arange(-0.5, 0.5, dx)
    x_2_stag = np.arange(-0.5, 0.5, dx)
    xx_1_stag, xx_2_stag = np.meshgrid(x_1_stag, x_2_stag)

    x_1_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    x_2_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    xx_1_c, xx_2_c = np.meshgrid(x_1_c, x_2_c)
    
    if PR_INT:
        mesh3 = np.vstack((xx_1_c[np.newaxis],xx_2_c[np.newaxis],np.zeros((N,N))[np.newaxis]))

    mu_0 = 0.25
    dt = 0.2*dx
    clock_max = np.int(np.round(20.0/dt,0))
    #clock_max = 1
    print 'grid size', N
    print 'total timesteps', clock_max

    #set up Lagrangian parameters
    #i_lag = np.arange(0, b_q)
    #ip_s_lag = np.remainder(i_lag+N_r, b_q)
    #im_s_lag = np.remainder(i_lag-N_r, b_q)
    k_teth = 0.5
    ks_teth = 10.0
    kb_teth = 0.01
    body_f = 0.5 
    #body_f = 0.0 
    rho_l = rho_0
    mu_l = mu_0

    #set up Lagrangian parameters
    N_r = 2*N
    dr = 1./N_r
    #X_1 = (a_lag+g_lag*(rr_lag-0.5))*np.cos(2*pi*ss_lag)
    #X_2 = (b_lag+g_lag*(rr_lag-0.5))*np.sin(2*pi*ss_lag)
    X_1_teth = np.arange(-.5, .5, dr)
    X_2_teth = .2*np.ones((N_r))
    X_1_teth = np.hstack((X_1_teth, np.arange(-.5, .5, dr))) 
    X_2_teth = np.hstack((X_2_teth, -.2*np.ones((N_r)))) 
    X_1_c = np.arange(-0.1+dx/4, 0.1+dx/4, dx/2)
    X_2_c = np.arange(-0.2+dx/4, 0.2+dx/4, dx/2)
    X_2_cu = np.arange(-0.1+dx/4, 0.1+dx/4, dx/2) 
    XX_1_c, XX_2_c = np.meshgrid(X_1_c, X_2_cu)
    N_s = XX_1_c.flatten(1).shape[0]
    X_1 = np.hstack((X_1_teth, .001*np.random.rand(N_s)+XX_1_c.T.flatten(1)))
    X_2 = np.hstack((X_2_teth, .001*np.random.rand(N_s)+XX_2_c.T.flatten(1)))
    b_q = X_1.shape[0]
    print b_q
    X_1 = X_1.reshape(b_q,1).copy()
    X_2 = X_2.reshape(b_q,1).copy()
    #X_1 = X_1_teth.reshape(b_q,1).copy()
    #X_2 = X_2_teth.reshape(b_q,1).copy()
    X_1_teth = X_1_teth.reshape(2*N_r,1).copy()
    X_2_teth = X_2_teth.reshape(2*N_r,1).copy()
    #Mu = 0.0*mu_0*(1.+np.sin(2*pi*rr_lag-pi/2))
    #Mu = 0.03*mu_0*np.ones((b_q,1))
    #Rho = 0.0*rho_0*(1.+np.sin(2*pi*rr_lag-pi/2))
    #Rho = 0.03*rho_0*np.ones((b_q,1))
    #Mu = np.zeros((b_q, 1)) 
    Mu = np.hstack((np.zeros((2*N_r)), mu_l*np.ones((N_s))))  
    Rho = np.hstack((np.zeros((2*N_r)), rho_l*np.ones((N_s))))  
    #Rho = np.zeros((b_q, 1)) 
    dtheta = dx**2/4.
    i_teth = np.arange(0, 2*N_r)
    i_t = np.arange(0, N_r)
    ip_t = np.hstack((np.arange(1,N_r),0))
    im_t = np.hstack((N_r-1,np.arange(0,N_r-1)))
    i_t = np.hstack((i_t, i_t+N_r))
    ip_t = np.hstack((ip_t, ip_t+N_r))
    im_t = np.hstack((im_t, im_t+N_r))
    body = (np.abs(xx_2_c)<0.2)*np.ones((N,N))
    N_r = 2*N_r

    #mu_bar = dtheta*np.max(Mu)
    #rho_bar = dtheta*np.max(Rho)
    mu_bar = 0.125 
    rho_bar = 0.5 
    #mu_bar = 0 
    #rho_bar = 0 
    print 'mu_bar is', mu_bar
    print 'rho_bar is', rho_bar
    mu_0 = mu_0+mu_bar
    rho_0 = rho_0+rho_bar

    firstn1 = -1*np.ones((N,N))
    nextn1 = -1*np.ones((b_q,1))
    firstn2 = -1*np.ones((N,N))
    nextn2 = -1*np.ones((b_q,1))

    start_time = time.time()
    for k in range(0,b_q):
        s_1_stag = (X_1[k]+0.5)/dx
        s_2_stag = (X_2[k]+0.5)/dx        
        i_1_stag = (np.floor(s_1_stag)).astype(np.uint32)
        i_2_stag = (np.floor(s_2_stag)).astype(np.uint32)
        s_1_c = (X_1[k]+0.5-dx/2)/dx;
        s_2_c = (X_2[k]+0.5-dx/2)/dx;        
        i_1_c = (np.remainder(np.floor(s_1_c),N)).astype(np.uint32);
        i_2_c = (np.remainder(np.floor(s_2_c),N)).astype(np.uint32);
        if (firstn1[i_2_c,i_1_stag]==-1):
            firstn1[i_2_c,i_1_stag]=k 
        else:
            nextn1[k]=firstn1[i_2_c,i_1_stag]
            firstn1[i_2_c,i_1_stag]=k
        if (firstn2[i_2_stag,i_1_c]==-1):
            firstn2[i_2_stag,i_1_c]=k 
        else:
            nextn2[k]=firstn2[i_2_stag,i_1_c]
            firstn2[i_2_stag,i_1_c]=k
    print 'time spent binning points is', time.time()-start_time, 'seconds'

    L_tf = calc_transf(dx)
    i_1, i_2, ip_1, ip_2, im_1, im_2 = calc_inds(N)

    #u_1_soln = 2*np.cos(2*pi*xx_1_stag)*np.sin(4*pi*xx_2_c)
    #u_2_soln = -np.sin(2*pi*xx_1_c)*np.cos(4*pi*xx_2_stag)
    u_1_soln = np.zeros((N,N)) #initial conditions
    u_2_soln = np.zeros((N,N))
    #u_1_soln = np.sin(2*pi*xx_2_c) #shear field
    #u_2_soln = np.zeros((N,N))

    #move data to gpu
    start_time = time.time()
    plan = Plan((N, N), queue=queue, dtype=np.complex128, fast_math=False)
    print 'time spent making plan is', time.time()-start_time, 'seconds'

    u_1_gpu = cl_array.to_device(queue, u_1_soln.astype(np.float64))
    u_2_gpu = cl_array.to_device(queue, u_2_soln.astype(np.float64))
    u_1_o_gpu = cl_array.to_device(queue, u_1_soln.astype(np.float64))
    u_2_o_gpu = cl_array.to_device(queue, u_2_soln.astype(np.float64))
    u_1_tf_gpu = cl_array.to_device(queue, u_1_soln.astype(np.complex128))
    u_2_tf_gpu = cl_array.to_device(queue, u_2_soln.astype(np.complex128))
    L_tf_gpu = cl_array.to_device(queue, L_tf.astype(np.complex128))
    firstn1_gpu = cl_array.to_device(queue, firstn1.astype(np.int32))
    nextn1_gpu = cl_array.to_device(queue, nextn1.astype(np.int32))
    firstn2_gpu = cl_array.to_device(queue, firstn2.astype(np.int32))
    nextn2_gpu = cl_array.to_device(queue, nextn2.astype(np.int32))
    f_1_gpu = cl_array.empty_like(u_1_gpu) 
    f_2_gpu = cl_array.empty_like(u_1_gpu) 
    f_1_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    f_2_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    mu_1_gpu = cl_array.empty_like(u_1_gpu) 
    mu_2_gpu = cl_array.empty_like(u_1_gpu) 
    if SAVE:
        rho_1_save_gpu = cl_array.empty_like(u_1_gpu) 
        rho_2_save_gpu = cl_array.empty_like(u_1_gpu) 
        mu_1_save_gpu = cl_array.empty_like(u_1_gpu) 
        mu_2_save_gpu = cl_array.empty_like(u_1_gpu) 
    if PR_INT:
        vort_gpu = cl_array.empty_like(u_1_gpu) 
    rho_1_gpu = cl_array.empty_like(u_1_gpu) 
    rho_2_gpu = cl_array.empty_like(u_1_gpu) 
    p_gpu = cl_array.empty_like(u_1_gpu) 
    p_last_gpu = cl_array.empty_like(u_1_gpu) 
    p_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_1_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    u_2_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    X_1_gpu = cl_array.to_device(queue, X_1.astype(np.float64))
    X_2_gpu = cl_array.to_device(queue, X_2.astype(np.float64))
    X_1_tilde_gpu = cl_array.to_device(queue, X_1.astype(np.float64))
    X_2_tilde_gpu = cl_array.to_device(queue, X_2.astype(np.float64))
    #rr_lag_gpu = cl_array.to_device(queue, rr_lag.astype(np.float64))
    #ip_s_lag_gpu = cl_array.to_device(queue, ip_s_lag.astype(np.uint32))
    #im_s_lag_gpu = cl_array.to_device(queue, im_s_lag.astype(np.uint32))
    F_1_gpu = cl_array.empty_like(X_1_gpu) 
    F_2_gpu = cl_array.empty_like(X_1_gpu) 
    Mu_gpu = cl_array.to_device(queue, Mu.astype(np.float64)) 
    Rho_gpu = cl_array.to_device(queue, Rho.astype(np.float64)) 
    U_1_gpu = cl_array.empty_like(X_1_gpu) 
    U_2_gpu = cl_array.empty_like(X_1_gpu) 
    X_1_teth_gpu = cl_array.to_device(queue, X_1_teth.astype(np.float64))
    X_2_teth_gpu = cl_array.to_device(queue, X_2_teth.astype(np.float64))
    i_teth_gpu = cl_array.to_device(queue, i_teth.astype(np.uint32))
    i_t_gpu = cl_array.to_device(queue, i_t.astype(np.uint32))
    ip_t_gpu = cl_array.to_device(queue, ip_t.astype(np.uint32))
    im_t_gpu = cl_array.to_device(queue, im_t.astype(np.uint32))
    body_gpu = cl_array.to_device(queue, body.astype(np.float64)) 

    #solve
    queue.finish()
    start_time = time.time()
    for t_step in range(1, clock_max+1):
        solve(t_step)
        if t_step % 50 == 0:
            print t_step
            queue.finish()
            if SAVE: 
                #print u_1_gpu.get().max()
                print 'u_1 max is', u_1_gpu.get().max()
                print 'u_2 max is', u_2_gpu.get().max()
                print 'rho_1 max is', rho_1_save_gpu.get().max()
                print 'rho_1 min is', rho_1_save_gpu.get().min()
                print 'mu_1 max is', mu_1_save_gpu.get().max()
                print 'mu_1 min is', mu_1_save_gpu.get().min()
                #print 'rho_1 max is', rho_1_save_gpu.get().max()
                #print 'rho_2 max is', rho_2_save_gpu.get().max()
        if PR_INT:
            if t_step % PR_INT == 0:
                from os.path import join as pjoin
                #vis_p = p_gpu.get().reshape(N,N,1)
                #file_name = "p-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("p", vis_p)])

                #vis_mu = mu_0+mu_1_save_gpu.get().reshape(N,N,1)
                #vis_u_1 = u_1_gpu.get().reshape(N,N,1)
                #file_name = "u_1-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("u_1", vis_u_1)])

                #vis_u_2 = u_2_gpu.get().reshape(N,N,1)
                #file_name = "u_2-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("u_2", vis_u_2)])

                vis_mu_2 = mu_2_save_gpu.get().reshape(N,N,1)
                file_name = "mu_2-%04d.vts"%t_step
                path_to_file = pjoin("vis_files", file_name) 
                write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                point_data=[("mu_2", vis_mu_2)])

                #vis_rho = rho_0+rho_1_gpu.get().reshape(N,N,1)
                #file_name = "rho-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("rho", vis_rho)])

                vis_vort = vort_gpu.get().reshape(N,N,1)
                file_name = "vort-%04d.vts"%t_step
                path_to_file = pjoin("vis_files", file_name) 
                write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                point_data=[("vort", vis_vort)])

                vis_X_1 = np.mod(X_1_gpu.get()+0.5, 1)-0.5
                vis_X_2 = np.mod(X_2_gpu.get()+0.5, 1)-0.5
                vis_X_3 = np.zeros((vis_X_1.shape))
                vis_X = np.hstack((vis_X_1,vis_X_2,vis_X_3)).copy()
                data = [
                ("F_1", F_1_gpu.get().reshape(b_q).copy()),
                ("F_2", F_2_gpu.get().reshape(b_q).copy()),
                ]

                file_name = "X-%04d.vtu"%t_step
                compressor = None

                grid = UnstructuredGrid(
                     (b_q, DataArray("points", vis_X, vector_format=VF_LIST_OF_VECTORS)),
                     cells=np.arange(b_q, dtype=np.uint32),
                     cell_types=np.asarray([VTK_VERTEX] * b_q, dtype=np.uint8))

                for name, field in data:
                   grid.add_pointdata(DataArray(name, field,
                   vector_format=VF_LIST_OF_COMPONENTS))

                from os.path import exists
                path_to_file = pjoin("vis_files", file_name)
                if exists(path_to_file):
                   raise RuntimeError("output file '%s' already exists"
                   % path_to_file)

                outf = open(path_to_file, "w")
                AppendedDataXMLGenerator(compressor)(grid).write(outf)
                outf.close()

    queue.finish()
    print 'time spent is', time.time()-start_time, 'seconds'

    u_1_soln = u_1_gpu.get()
    u_2_soln = u_2_gpu.get()
    X_1_soln = X_1_gpu.get()
    X_2_soln = X_2_gpu.get()
    p_soln = 0.5*(p_gpu.get()+p_last_gpu.get())
    #print u_1_soln
    #print u_2_soln
    #calculate error
    if ii != start_size: 
        errinf1, errinf2, err11, err12, err21, err22, errinfX, err1X, err2X, errinfp, err1p, err2p = err(u_1_soln, u_2_soln, u_1_o, u_2_o, X_1_soln, X_2_soln, X_1_o, X_2_o, p_soln, p_o)
        #errinf1, errinf2, err21, err22 = err(u_1_soln, u_2_soln, u_1_o, u_2_o)
        print 'u Errinf1 is', errinf1, '\n u Errinf2 is', errinf2
        print 'u Err11 is', err11, '\n u Err12 is', err12
        print 'u Err21 is', err21, '\n u Err22 is', err22
        print 'X Errinf is', errinfX
        print 'X Err1 is', err1X
        print 'X Err2 is', err2X
        print 'p Errinf is', errinfp
        print 'p Err1 is', err1p
        print 'p Err2 is', err2p
        if ii > start_size+1:
            print 'rinf for u1', -np.log(errinf1/errinf1_o)/np.log(2)
            print 'rinf for u2', -np.log(errinf2/errinf2_o)/np.log(2)
            print 'r1 for u1', -np.log(err11/err11_o)/np.log(2)
            print 'r1 for u2', -np.log(err12/err12_o)/np.log(2)
            print 'r2 for u1', -np.log(err21/err21_o)/np.log(2)
            print 'r2 for u2', -np.log(err22/err22_o)/np.log(2)
            print 'rinf for X', -np.log(errinfX/errinfX_o)/np.log(2)
            print 'r1 for X', -np.log(err1X/err1X_o)/np.log(2)
            print 'r2 for X', -np.log(err2X/err2X_o)/np.log(2)
            print 'rinf for p', -np.log(errinfp/errinfp_o)/np.log(2)
            print 'r1 for p', -np.log(err1p/err1p_o)/np.log(2)
            print 'r2 for p', -np.log(err2p/err2p_o)/np.log(2)
        errinf1_o = errinf1
        errinf2_o = errinf2
        err11_o = err11
        err12_o = err12
        err21_o = err21
        err22_o = err22
        errinfX_o = errinfX
        err1X_o = err1X
        err2X_o = err2X
        errinfp_o = errinfp
        err1p_o = err1p
        err2p_o = err2p
    u_1_o = u_1_soln
    u_2_o = u_2_soln
    X_1_o = X_1_soln
    X_2_o = X_2_soln
    p_o = p_gpu.get()
